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Abstract 

This  document  constitutes  the  final  project  report  for  Contract  #  FA9550-04-1- 
0187  titled  Theory  and  Applications  of  Computational  Time-reversal  Imaging.  The 
report  summarizes  the  theoretical  development  and  implementation  of  time-reversal 
based  imaging  and  target  detection  algorithms  for  locating  targets  from  multistatic 
data  collected  from  unstructured  phased  antenna  arrays.  The  report  summarizes  the 
main  theoretical  results  obtained  in  the  program  and  includes  both  computer  simulated 
examples  as  well  as  results  from  experimental  data  collected  by  a  research  team  from 
Carnegie  Mellon  University  illustrating  the  use  of  the  algorithms  developed  in  the 
project.  The  final  section  of  the  report  outlines  goals  for  follow-on  research  in  these 
general  areas. 
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1  Project  Summary 

The  NU  team  consisted  of  myself  as  PI  as  well  as  Professors  Hanoch  Lev-Ari  and  Edwin 
Marengo  (NU),  Professor  Arye  Nehorai  (UI)  and  Drs.  Sean  Lehman  and  David  Chambers 
(LLNL).  During  this  project  the  NU  team  was  directed  by  DARPA  to  concentrate  most  of 
its  effort  on  the  problem  of  detecting  a  point  target  embedded  in  heavy  clutter.  They  have 
modeled  this  problem  exactly  using  the  well-known  Foldy  Lax  equations  [1]  that  incorporate 
all  multiple  scattering  between  the  targets,  clutter  and  background  and  have  proven  [2]  that 
the  (non-linear)  detection  problem  can  be  decomposed  into  the  two  separate  problems  of 

1.  Estimating  the  target/clutter  locations  using  either  classical  time-reversal  imaging 
(DORT)  [11]  or  time  reversal  MUSIC  [3,  4], 

2.  Estimating  the  target/clutter  scattering  amplitudes  using  a  non-linear  iterative  algo¬ 
rithm  developed  in  the  study  [5] . 

It  is  important  to  emphasize  that  even  in  the  presence  of  multiple  scattering  the  above  two 
problems  are  independent;  i.e.,  the  target/clutter  locations  can  be  estimated  using  either 
time  reversal  DORT  or  MUSIC  without  knowing  their  scattering  amplitudes.  Once  the 
target  and  clutter  locations  are  known  their  scattering  amplitudes  are  determined  using  a 
non-linear  iterative  algorithm  that  was  also  developed  within  the  course  of  the  project  [5]. 
The  final  step  in  the  detection  problem  is  to  use  the  computed  target/clutter  amplitudes  in 
a  maximum  likelihood  based  detection  algorithm.  This  time  reversal  based  detection  scheme 
has  been  tested  and  evaluated  in  extensive  Monte  Carlo  computer  simulations  against  other 
(linear)  time  reversal  based  detection  schemes  and  against  benchmark  matched  filter  based 
schemes  [6] .  In  all  cases  tested  the  time  reversal  approaches  always  outperform  the  matched 
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filter  based  methods  and  the  non-linear  time  reversal  based  scheme  developed  in  the  project 
always  came  out  best. 

The  target/clutter  location  estimation  step  is  of  key  importance  in  the  detection  problem 
and  the  NU  team  has  spent  a  good  deal  of  effort  on  developing  and  testing  alternative  time 
reversal  based  imaging  algorithms  for  performing  this  step.  An  important  breakthrough  in 
this  investigation  has  recently  been  achieved  when  it  was  shown  that  time  reversal  MUSIC 
is  only  one  of  a  large  class  of  non-linear  algorithms  that  allow  point  target/clutter  locations 
to  be  estimated  from  the  singular  vectors  of  the  SVD  of  the  multistatic  data  matrix  even  in 
the  presence  of  intense  multiple  scattering  between  targets  and  clutter.  These  generalized 
MUSIC  algorithms  are  now  being  tested  and  evaluated  on  experimental  data  generated  by  a 
reserch  team  from  Carnegie  Mellon  University  (CMU)  and  a  journal  paper  jointly  authored 
by  NU  and  CMU  team  members  on  this  work  is  in  preparation  [7]. 

Besides  developing,  testing  and  evaluating  time  reversal  based  detection  algorithms  the 
NU  team  has  also  spent  a  good  deal  of  effort  developing  time  reversal  based  inverse  scat¬ 
tering  algorithms  that  yield  quantitative  estimates  of  distributed  scatterers  from  multistatic 
scattered  field  data  [8,  9,  10].  Thus,  where  time  reversal  DORT  [11]  and  MUSIC  [3,  4]  and 
similar  algorithms  allow  sets  of  isolated  point  scatterers  to  be  imaged,  the  inverse  scattering 
algorithms  allow  extended  (3D)  targets  as  well  as  the  distributed  background  medium  to 
be  quantitatively  determined.  The  estimation  of  backgrounds  from  scattered  field  data  is 
extremely  important  within  the  context  of  the  project  since  knowledge  of  the  background 
allows  the  background  Green  function  to  be  readily  computed  and  used  in  the  time  rever¬ 
sal  based  detection  algorithms  developed  in  the  study.  It  is  important  to  mention  that  the 
inverse  scattering  algorithms  developed  in  the  study  are  based  on  the  SVD  of  the  mapping 
from  the  medium  to  the  multistatic  data  and,  hence,  are  extremely  fast  and  efficient  and 
can  thus  be  incorporated  into  near  real  time  applications. 

An  important  part  of  the  NU  research  was  the  development  of  more  conventional  ap¬ 
proaches  to  the  target  detection  problem  [12,  13].  In  particular,  a  number  of  maximum 
likelihood  (ML)  based  detection  schemes  using  standard  signal  processing  techniques  based 
on  the  time  reversal  scattering  model  for  the  multistatic  data  matrix  were  developed  and 
evaluated  in  a  Monte  Carlo  simulation  study.  One  of  these  schemes,  based  on  a  single  target 
scattering  model,  was  employed  as  the  benchmark  algorithm  for  evaluation  of  the  time  re¬ 
versal  based  algorithms  developed  in  the  program.  This  work  was  important  since  it  showed 
that  the  time  reversal  based  detection  algorithms  could  be  incorporated  into  the  mathemat¬ 
ical  framework  of  conventional  target  detection  theory  and  also  because  it  laid  the  ground 
work  of  the  extension  of  the  time  reversal  based  schemes  to  the  stochastic  framework  that 
will  be  employed  in  Phase  II  of  the  project  (see  below). 

2  Time  reversed  DORT  and  MUSIC  Images 

Much  of  the  NU  team  effort  was  directed  at  imaging  which  is  a  first  and  crucial  step  in  the 
target  detection  process.  In  imaging  the  acquired  multistatic  data  is  processed  in  such  a 
way  as  to  generate  images  of  a  set  of  targets  from  the  scattered  field  data.  We  used  three 
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schemes  for  generating  these  images: 

1.  Basic  time  reversal 

2.  The  DORT  method 

3.  Time  reversal  MUSIC 

These  three  methods  are  described  extensively  in  the  literature  cited  at  the  end  of  the 
report.  Images  from  both  simulated  as  well  as  experimental  data  provided  by  Carnegie 
Mellon  University  (CMU)  were  employed  in  the  study.  In  fact,  the  NU  team  were  the  first  to 
successfully  image  the  first  CMU  data  set  which  was  reported  in  a  preliminary  report  titled 
“Basic  Time  Reversal  Imaging  of  CMU  Data  Set,”  and  which  was  issued  early  in  the  first 
year  of  the  project.  In  the  remainder  of  this  section  we  first  summarize  the  major  results  of 
that  report  and  then  present  results  obtained  from  subsequent  CMU  data  sets  obtained  near 
the  end  of  the  project.  All  coding  was  done  in  Matlab  and  the  codes  with  accompanying 
written  reports  were  placed  on  the  DARPA  web  site  during  the  course  of  the  project. 

2.1  Early  Results  from  CMU  experimental  data 

Carnegie  Mellon  University  (CMU)  provided  experimental  test  data  on  which  the  DARPA 
teams  were  required  to  apply  their  target  imaging  and  detection  algorithms.  The  experi¬ 
mental  configuration  is  shown  in  fig.  1  and  consists  of  a  number  of  cylindrical  targets  whose 
radii  were  on  the  order  of  the  EM  wavelength  used  in  the  experiments.  Multistatic  data  were 
acquired  and  made  available  to  the  various  team  members  for  test  and  evaluation  of  their 
algorithms.  In  the  report  cited  above  we  used  “basic  time  reversal  imaging”  to  generate  im¬ 
ages  of  target  sets  from  the  first  CMU  data  set  acquired  in  the  project.  We  first  summarize 
how  basic  time  reversal  imaging  is  performed  and  then  present  the  imaging  results  obtained 
in  that  report. 

2.1.1  Basic  Time  Reversal  Imaging 

The  basic  method  of  time  reversal  imaging  is  simply  to  time  reverse  the  scattered  field 
measurements  from  all  array  sensor  locations.  At  any  given  frequency  this  corresponds 
to  the  process  of  field  back  propagation  in  the  background  medium  and  is  accomplished  by 
conjugating  the  measured  field  amplitudes  over  the  sensor  locations  and  propagating  the 
resulting  amplitudes  into  the  region  containing  the  scatterers  using  the  background  Green 
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function.  Mathematically,  this  generates  the  time  reversed  (back  propagated)  field 


where 
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are  the  transmit  and  receive  array  point  spread  functions  (PSF’s)  ,  respectively1.  At  any 
given  frequency  the  image  field  ijj  formed  from  basic  time  reversal  imaging  from  a  set  of 
point  scatterers  is  seen  to  be  equal  to  a  sum  of  the  product  of  the  array  PSF’s  weighted  by 
the  scatterers’  reflection  coefficients  all  centered  at  the  scatterer  locations. 


2.1.2  Time  Domain  Field 

The  product  of  the  transmit  and  receive  array  PSF’s  //^(x,  xm,  w)//b(x,  xm,  uj)  will  peak  at 
the  scatterer  locations  x  =  xm  so  that  the  peaks  in  the  magnitude  of  the  time  reversed  field 
ip  give  an  estimate  of  the  locations  of  the  various  scatterers;  i.e.,  will  peak  at  the  various 
scatterer  locations  xm,  m  =  1, 2,  •  •  •  ,  M.  The  resolution  of  this  basic  scheme  for  estimating 
target  location  can  be  improved,  especially,  in  backscatter  measurements  such  as  is  the  case 
for  the  CMU  data,  by  integrating  the  field  amplitude  over  the  frequency  band  of  the  data. 
In  particular, 

1  f°° 

Ha,b  (x,  xm,  t)  =  —  J  do)  Ha,b  (x,  xm,  u)e~lujt 

will  achieve  its  peak  value  at  t  =  0  so  that  the  time  domain  time  reversed  field  evaluated  at 
t  =  0  will  yield  optimum  resolution  of  the  various  scatterer  locations.  This  quantity  is  given 

by 

i  r°°  l  r°°  jE ' 

®(x)  =  —  /  dunp(x,u))  =  —  /  (L)J2J2Klk(UJ)Go(f3j,yi,u))G0(ak,x,uj). 

Z7f  J -oo  Z7F  j= 1  fc=l 

1Note  that  Ha  only  involves  the  background  Green  function  while  Hg  involves  both  the  background  as 
well  as  the  full  Green  functions.  It  can  be  shown  that  the  peaks  of  both  G0  and  G  are,  however,  always 
occur  at  the  point  x  =  xm;  i.e.,  at  the  target  location. 
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Up  Na 

*w  =  EEE  KlkMGoiPj,*,  u>n)Go(ack,  x,  un) 

Wn  j=  1  k=  1 

2.2  Processed  Image  Data  from  first  CMU  data  set 

We  processed  all  of  the  first  CMU  data  set  using  both  the  basic  time  reversal  imaging 
scheme  described  above  as  well  as  an  SVD  scheme  described  in  the  early  report  and  useful 
for  well  resolved  scatterers.  We  present  in  figures  2-5  the  results  obtained  using  the  basic 
time  reversal  imaging  scheme.  The  images  were  generated  using  the  code  cmul.m  which  we 
placed  on  the  DARPA  web  site. 

2.3  Processed  Images  from  later  CMU  data  sets 

We  show  in  fig.  6  a  schematic  diagram  of  the  CMU  experimental  configuration  used  in  a 
number  of  the  later  tests.  In  data  set  #  11  a  total  of  17  dielectric  rods  were  used  as  targets 
while  in  data  set  #  12  a  conducting  rod  was  added  to  the  configuration.  We  imaged  first 
the  dielectric  rods,  then  the  dielectric  rods  plus  conducting  rod  then  finally  processed  the 
difference  data  to  generate  a  final  image  of  the  conducting  rod  alone.  The  imaging  results 
for  standard  time  reversal  imaging  are  shown  in  fig.  7  and  those  obtained  using  time  reversal 
MUSIC  are  shown  in  fig.  8.  The  two  are  compared  in  fig.  9. 
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Time  Reversed  Field 


Contour  Plot  of  Time  Reversed  Field 


Figure  3:  Time  reversed  field  from  eight  scatterers  processed  at  ten  frequencies  equally 
spaced  between  4.5  and  5.5  MHz.  The  number  of  the  scatterer  is  illustrated  in  both  figures 
and  the  x,y  axes  are  in  centimeters.  Note  that  due  to  the  quirks  of  Matlab  that  the  contour 
plot  is  drawn  on  an  inverted  vertical  axis. 
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Figure  4:  Time  reversed  field  from  ten  scatterers  processed  at  ten  frequencies  equally  spaced 
between  4.5  and  5.5  MHz.  The  number  of  the  scatterer  is  illustrated  in  both  figures  and  the 
x,y  axes  are  in  centimeters.  Note  that  due  to  the  quirks  of  Matlab  that  the  contour  plot  is 
drawn  on  an  inverted  vertical  axis. 


Time  Re\ersed  Field 


Contour  Plot  of  Time  Reversed  Field 


Figure  5:  Time  reversed  field  from  fourteen  scatterers  processed  at  twenty  frequencies  equally 
spaced  between  4.5  and  5.5  MHz.  The  number  of  the  scatterer  is  illustrated  in  both  figures 
and  the  x,y  axes  are  in  centimeters.  Note  that  due  to  the  quirks  of  matlab  that  the  contour 
plot  is  drawn  on  an  inverted  vertical  axis. 
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CMU  3rd  Data  Sets  11  &12 


17  dielectric  rods 


17  dielectric  rods  +  conducting  rod 


Processed  with  standard  Time  Reversal  Imaging  as  well  as  with  TR  MUSIC 

April  1 5, 2005  NU  Team  Darpa  Report 


Figure  6:  Scatterer  and  Antenna  configurations  for  CMU  experimental  data  for  CMU  Set 
#3  Setups  11  and  12. 
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CMU  3rd  Data  Sets  “Setupll”  and  ''Setup12” 
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Total  Time  Reversed  Image  from  difference  data 

Processed  with  7  frequencies  unifoimiy  spaced  between  4.84  and  5.14  Ghz 
April  1 5, 2005  NU  Team  Darpa  Report 


Figure  7:  Total  time  reversed  images  from  Experimental  Set  #3  Setups  11  and  12. 
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CMU  3rd  Data  Sets  “Setupl  1”  and  “Setup  12” 


Time  Reversal  MUSIC  Image  from  deleclric  rods  plus  conductor 


-100  -80  -60  -40  -20  0  20  40  60  60  100 

Time  Reversal  MUSIC  Image  from  difference  data 

Processed  with  7  frequencies  uniformly  spaced  between  4.84  and  5.14  Ghz 

April  15,2005  NU  Team  Darpa  Report 


Figure  8:  Time  reversal  MUSIC  images  from  Experimental  Set  #3  Setups  11  and  12. 
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CMU  3rd  Data  Sets  “Setupl  1”  and  “Setup  12” 
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Time  Reversal  Image  fromdelectric  rods  plus  conductor 
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Processed  with  7  frequencies  uniformly  spaced  between  4.84  and  5.14  Ghz 

April  1 5, 2005  NU  Team  Darpa  Report 


Figure  9:  Comparison  of  Time  reversal  and  Time  reversal  MUSIC  images  from  Experimental 
Set  #3  Setups  11  and  12. 
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3  DARPA  Test  Simulations 

As  mentioned  earlier,  one  of  the  main  objectives  of  this  program  was  the  detection  of  single 
or  multiple  point  targets  embedded  in  heavy  clutter.  The  NU  team  employed  an  exact 
scattering  model  for  this  problem  and  first  evaluated  this  model  in  computer  simulations. 
We  present  in  this  section  the  results  of  the  earliest  such  study  which  was  a  Monte  Carlo  test 
simulation  consisting  of  the  detection  of  a  single  point  target  embedded  in  the  presence  of  a 
number  of  point  clutter  targets.  The  simulation  tests  the  performance  of  two  time  reversal 
based  detection  algorithms;  one  based  on  the  use  of  the  distorted  wave  Born  approximation 
(DWBA)  for  the  K  matrix  and  a  second  using  a  full  multiple  scattering  based  model  for  the 
K  matrix  [2,  5].  The  details  of  the  simulation  code  are  reviewed  and  a  few  examples  of  the 
ROC  and  other  performance  curves  resulting  from  the  code  are  presented. 

3.1  Synthetic  Data  Generation 

We  wish  to  compute  a  single  realization  of  the  measurement  of  the  K  matrix  for  a  system  of 
one  dominant  point  target  in  the  presence  of  Mc  <  N  —  1  clutter  targets  within  a  background 
whose  Green  function  G0  is  known.  We,  of  course,  want  to  include  all  multiple  scattering 
between  the  scatterers  in  the  computation.  The  locations  xm,  m  —  0, 1,  •  •  •  ,  Mc  of  the  target 
and  clutter  are  random  variables  uniformly  distributed  within  a  square  with  edge  length  Lx 
which  is  a  free  parameter  in  the  simulation  code.  The  scattering  amplitude  of  the  target 
To  is  randomly  chosen  to  be  either  0  or  unity  (equally  probable)  and  the  clutter  scattering 
amplitudes  rm,m  =  1,2,  •••  ,MC  are  taken  to  be  zero  mean  Gaussian  having  a  standard 
deviation  am  that  is  also  a  free  parameter.  The  simulation  code  works  at  a  single  frequency 
u>  which  will  be  assumed  to  be  fixed  throughout  the  following  discussion.  The  frequency  is 
fixed  in  the  code  by  the  use  of  a  background  wavelength  A  which  is  set  to  unity.  All  lengths 
in  the  code  are  then  given  in  wavelength  units. 

For  any  given  realization  of  the  scattering  center  locations  {xm}  and  scattering  and  target 
strengths  {rm}  the  code  generates  synthetic  data  using  the  Foldy  Lax  equations  [1] 

G(xm,ak,to)  =  Go(xm,ak,u)  +  rm'(w)G0(xm,xm',w)G(xm',afc,w)  (3a) 

m'^m 

where  the  antenna  element  locations  ak,  k  =  1,2,-  ••  ,  N  are  also  free  parameters  in  the 
simulation  code.  The  Foldy  Lax  equations  are  solved  iteratively  and  the  resulting  full  Green 
function  G(xm,  ak.  u>)  is  then  used  to  compute  the  (noise  free)  multistatic  data  matrix  ac¬ 
cording  to  its  definition 

Mc 

KM  =  J2Tm  (w)Go  (cKj ,  Xjjj,  Ld)G{xm,  OLk ,  ui)  (3b) 

m=0 

The  final  step  in  the  data  generation  process  is  the  addition  of  additive  noise  to  the 
computed  multistatic  data  matrix  according  to  the  equation 

Kj>k(uj)  =  Kj.fcH  +  Max  |Rj,fc|A/j,fc(0,  aw)el<i>i'k  (4) 
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where  A/j,fc(0,  aw)  is  a  random  matrix  whose  elements  are  uncorrelated  and  zero  mean  Gaus¬ 
sian  with  standard  deviation  equal  to  aw  and  where  4>jtk  is  an  uncorrelated  random  matrix 
whose  elements  are  uniformly  distributed  between  0  and  2tt.  The  noisy  K  matrix  K  is  then 
input  into  the  time  reversal  based  target  detection  algorithms  described  below. 

3.2  Time  Reversal  MUSIC 

Both  of  the  time  reversal  based  algorithms  described  below  are  assumed  to  have  a  first 
step  that  consists  of  estimation  of  all  the  scattering  centers  {xm}  (this  includes  target  plus 
clutter).  We  showed  previously  that  as  long  as  the  number  of  scattering  centers  M  —  Mc  + 1 
is  less  than  the  number  of  antenna  elements  N  that  this  can  be  accomplished  with  time 
reversal  MUSIC  using  the  background  Green  function  and  the  singular  vectors  vp,  ap  =  0 
where  Kvp  =  apup  [3,  4,  15].  We  assume  here  that  this  step  is  performed  exactly  by  both  of 
the  algorithms  discussed  below. 

3.3  Scattering  Strength  Estimation 

Once  the  scattering  centers  are  estimated  (assumed  known  in  the  current  Monte  Carlo  sim¬ 
ulation  code)  the  scattering  strengths  rm,m  =  0,1,- ••  ,MC  are  estimated.  We  used  two 
procedures  for  doing  this:  one  based  on  the  DWBA  model  for  the  multistatic  data  matrix 
and  a  second  based  on  the  full  multiple  scattering  model  for  this  quantity.  Thus,  we  are,  in 
fact,  using  two  time  reversal  algorithms:  one  based  on  the  DWBA  and  a  second  on  the  full 
multiple  scattering  model  using  the  Foldy  Lax  equations.  We  compared  both  of  these  algo¬ 
rithms  with  a  standard  “benchmark”  algorithm  in  our  work  with  the  CMU  data  presented 
later  in  the  report. 

3.3.1  DWBA  Based  Algorithm 

The  DWBA  algorithm  employs  the  DWBA  model  for  the  multistatic  data  matrix  for  a  system 
of  N  antennas  located  at  atj,  j  =  1, 2,  •  •  •  ,  N  and  M  =  Mc  +  1  point  scatterers  located  at 
xm,  m  =  0, 1, 2,  •  •  •  ,  Mc,  where,  as  above,  Mc  is  the  number  of  clutter  targets: 

Mc 

Kjk  =  ^rTn(w)Go(aj,xm,w)G0(xm,afc,w),  (5) 

m= 0 

where  we  have  used  the  superscript  “b”  to  denote  the  DWBA  (Born)  model  for  the  K  matrix. 
We  write  the  above  set  of  equations  in  the  form 

Mc 

Kj,k  =  Go(ajtXm,u)Am(k)  (6a) 

771=0 

where 

Am(k)  =  rm(w)G0(xra,at,w).  (6b) 
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Eqs(6a)  are  a  set  of  N2  equations  for  the  N  x  M  unknowns  Am(k),  m  =  0, 1,  •  •  •  ,  Mc, 
k  =  1, 2,  •  •  •  ,N.  In  the  code  the  Kb-  k  is  set  equal  to  the  computed  noisy  K  matrix  Kj <k  and 
a  least  squares  solution  of  Eqs.(6a)  is  obtained  for  each  value  of  k  =  1, 2,  •,  N.  The  final  step 
is  the  estimation  of  the  scattering  coefficients  which,  in  the  DWBA  model,  are  computed 
using  the  algorithm 


£ 


dm(fc) 

G0(-xm,ak,uj) 


(7) 


which  is  simply  the  an  average  of  the  N  values  of  rm  given  according  to  Eq.(6b). 


3.3.2  Multiple  Scattering  Based  Algorithm 

The  multiple  scattering  based  time  reversal  detection  algorithm  uses  the  exact  multiple 
scattering  model  Eq.(3b)  for  the  multistatic  data  matrix.  This  algorithm  differs  from  the 
DWBA  algorithm  described  above  only  in  the  definition  of  the  quantities  Am(k )  and  in  the 
computation  of  the  scattering  strengths  from  the  least  squares  solution  of  the  set  Eqs.(6a) 
(with,  of  course,  Kbk  set  equal  to  the  noisy  measured  K  matrix  Kj,k).  In  particular,  in  place 
of  Eq.(6b),  the  amplitudes  Am(k)  in  the  multiple  scattering  model  are  defined  by 


J4m(k)  —  Tm(cj)G(xm,  Qlfc,  w). 


(8a) 


where  G  is  the  full  Green  function  and  must  satisfy  the  Foldy  Lax  equations  Eqs.(3a).  In 
terms  of  Am(k )  these  equations  can  be  written  in  the  form 

G(xm,c>:fc,cj)  =  G0(xm,afc,cu)  -f  Go(xm,xm/,w)Am/(fc).  (8b) 

m'^m 


The  multiple  scattering  based  algorithm  thus  first  uses  a  least  squares  solution  of  the  set 
Eqs.(6a)  which  are  then  input  into  the  Foldy  Lax  equations  Eqs.(8b)  to  compute  an  estimate 
G  of  the  full  Green  function  G(xm,  ack,u>).  The  scattering  strengths  are  then  estimated  via 
the  equation 


T’m 


J_  Am(k) 

N  G(xm,  ak,u>)  ’ 


(9) 


which  is  seen  to  differ  form  the  DWBA  solution  Eq.(7)  only  in  the  replacement  of  the 
background  Green  function  Go  with  the  (estimated)  full  Green  function  G. 


3.4  Detection  Algorithms 

The  first  processing  step  for  both  algorithms  is  target  location  estimation  via  MUSIC.  This 
is  a  non-linear  step  which  makes  the  remaining  estimation  step  of  the  scattering  strengths 
itself  non-linear  so  that  it  is  not  possible  to  determine  an  “optimum”  detection  algorithm. 
However,  if  we  make  the  simplifying  assumption  that  this  first  step  is  exact  (which  is  a 
reasonable  assumption  since  MUSIC  does  perform  very  well)  the  estimation  of  the  scattering 


18 


i 


A 


strengths  for  the  DWBA  algorithm  is  linear  and  an  optimum  algorithm  can  then  be  devised2. 
However,  the  same  is  not  true  for  the  multiple  scattering  based  algorithm  since  the  Foldy 
Lax  equations  are  non-linear  with  the  result  that  the  estimation  of  the  full  Green  function 
and  target  scattering  strengths  is  a  non-linear  process.  However,  for  the  sake  of  simplicity  we 
decided  to  use  a  simple  threshold  based  algorithm  for  the  detection  step  in  both  algorithms. 
Such  an  algorithm  would  only  be  appropriate  for  a  linear  model  with  AWGN  but  is  used 
here  for  lack  of  anything  better. 

3.4.1  Threshold  Test 

In  any  given  iteration  (realization  of  random  parameters)  the  Monte  Carlo  code  generates 
two  solution  vectors  f  =  [f0, n,---  ,tmc]  of  the  target  scattering  strengths  (one  using  the 
DWBA  algorithm  and  one  using  the  multiple  scattering  based  algorithm  described  above). 
The  actual  target  strength  is  To  and  the  remaining  estimates  are  of  the  clutter  amplitudes 
Tm,m  =  1,2,  •  •  •  ,  Mc.  Both  the  DWBA  as  well  as  the  multiple  scattering  based  algorithm 
make  a  decision  by  comparing  the  maximum  value  of  the  real  part  of  the  estimated  solution 
vectors  J  =  Max  9?(f)with  a  threshold  T  and  deciding  that  the  target  is  present  if  J  >  T 
and  not  present  if  J  <  T.  The  probabilities  of  detection  Pd  and  false  alarm  Pf  are  then 
computed  from  these  results  in  a  standard  fashion  by  keeping  track  of  correct  and  incorrect 
decisions.  Plots  of  Pd  and  Pf  as  a  function  of  the  threshold  value  T  are  then  generated  as  is 
a  plot  of  the  Receiver  Operation  Curve  (ROC)  which  is  simply  a  plot  of  Pd  versus  Pf. 

3.5  Examples 

To  illustrate  the  Monte  Carlo  program  we  show  a  few  examples  that  compare  the  performance 
of  the  DWBA  based  detection  algorithm  with  the  multiple  scattering  based  algorithm.  In 
these  and  later  studies  the  figures  are  presented  in  the  last  section  of  the  report.  In  all  of 
these  examples  a  linear  array  of  N  =  6  unequally  spaced  antenna  elements  were  employed 
and  Mc  =  4  clutter  targets  and  one  single  real  target  were  used.  The  target  and  clutter 
locations  were  uniformly  distributed  within  a  square  of  length  Lx  which  is  labeled  on  the 
figures  presented  in  the  last  section  of  the  report.  The  center  of  the  target/clutter  support 
square  was  32  wavelengths  from  the  antenna  array. 

In  Fig.  10  we  show  curves  of  the  probabilities  of  detection  Pd  and  false  alarm  Pf  plotted 
versus  threshold  value  T  for  one  set  of  simulation  parameters.  In  these  and  the  following 
examples,  the  target  has  a  scattering  strength  r0  that  is  equally  probable  to  be  unity  (target 
present)  or  zero  (no  target)  and  the  clutter  targets  were  independent  Gaussian  random 
variables  having  zero  mean  and  standard  deviation  oT  that  is  labeled  on  the  figures.  Also 
labeled  on  the  figures  are  the  number  of  Monte  Carlo  runs,  the  noise  standard  deviation 
aw  (see  Eq.(4)),  and  the  length  of  the  side  of  the  support  square  for  targets  and  clutter  in 

2This  is  possible  since  both  the  additive  noise  as  well  as  the  clutter  amplitudes  are  Gaussian.  This 
algorithm  will,  of  course,  be  optimum  only  within  the  DWBA  model  and  does  not  account  for  multiple 
scattering  effects. 
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number  of  wavelengths.  The  ROC  curves  corresponding  to  the  results  shown  in  Fig.  10  are 
presented  in  Fig.  11. 

It  is  seen  from  Fig.  10  that  the  probability  of  detection  differs  significantly  for  the  two 
algorithms  but  that  the  false  alarm  probabilities  are  almost  identical.  This  same  trend  was 
found  in  all  simulations  that  were  performed  and  so  we  only  show  the  ROC  curves  for  the 
later  examples.  It  is  clear  from  both  Figs.  10  and  11  that  the  multiple  scattering  based 
algorithm  clearly  out  performs  the  DWBA  based  algorithm.  This  is  due  to  the  fact  that  the 
targets  and  clutter  are  closely  packed  (within  a  square  having  sides  of  length  Lx  =  3.2  A  and 
that  the  clutter  amplitudes  are  relatively  large  (ay  =  0.4). 

The  ROC  curves  for  the  same  set  of  parameters  as  used  in  the  first  example  but  with 
a  different  additive  noise  standard  deviation  are  shown  in  Fig.  12.  The  increase  in  additive 
antenna  noise  clearly  degrades  the  performance  of  both  algorithms  but  is  most  strongly  felt 
by  the  multiple  scattering  algorithm.  This  is  probably  due  to  error  introduced  into  the 
estimation  of  the  full  Green  function  from  the  set  of  Eqs.(8b)  due  to  error  in  the  coefficients 
Am(k).  Other  simulations  (not  shown)  indicate  that  this  trend  increases;  i.e.,  the  multiple 
scattering  based  algorithm  is  more  sensitive  to  additive  antenna  noise  than  is  the  DWBA 
based  algorithm. 

In  the  final  example  the  additive  noise  is  again  reduced  to  the  level  of  the  first  example 
but  the  standard  deviation  of  the  clutter  target  amplitudes  is  increased.  Again  the  multiple 
scattering  algorithm  out  performs  the  DWBA  based  algorithm  but  the  performance  of  both 
algorithms  degrade  relative  to  the  performance  in  the  first  example.  This  is  obviously  due 
to  false  estimates  of  target  presence  due  to  the  use  of  a  simple  threshold  based  detection 
algorithm;  i.e.,  the  maximum  value  of  the  estimated  r  vector  is  more  likely  due  to  a  high 
amplitude  clutter  target. 

4  Monte  Carlo  Simulation  from  CMU  Experimental 
Data 

Carnegie  Mellon  University  (CMU)  provided  experimental  test  data  on  which  the  DARPA 
teams  were  required  to  apply  their  target  detection  algorithms.  The  experimental  configu¬ 
ration  is  shown  in  fig.  1  and  consists  of  a  number  of  cylindrical  targets  whose  radii  were  on 
the  order  of  the  EM  wavelength  used  in  the  experiments.  Multistatic  data  were  acquired 
and  made  available  to  the  various  team  members  for  test  and  evaluation  of  their  algorithms. 

4.1  Signal  Model 

The  third  CMU  data  set  was  used  to  generate  ROC  curves  using  a  change  detection  based 
signal  model.  In  this  study  we  used  the  experimentally  obtained  multistatic  data  matrix 
KSj\c{id)  corresponding  to  a  given  set  of  clutter  targets  plus  a  single  copper  target  and  the 
experimentally  obtained  multistatic  data  matrix  K^k{uS)  corresponding  to  only  the  clutter 
targets  at  the  single  frequency  of  5  Ghz.  The  study  assumed  that  a  noisy  version  of  the 
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False  Alarm  Probability  P,  Probability  of  Detection  P 


Probability  of  Detection  versus  threshold 


Threshold  T 


Threshold  T 

Figure  10:  Plots  of  Probabilities  of  detection  and  false  alarm  for  the  DWBA  (dotted)  and 
multiple  scattering  based  (solid)  algorithms  versus  threshold  value. 
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Figure  11:  Receiver  Operating  Curves  for  the  DWBA  (dotted)  and  multiple  scattering  based 
(solid)  algorithms. 
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Figure  12:  Receiver  Operating  Curves  for  the  DWBA  (dotted)  and  multiple  scattering  based 
(solid)  algorithms. 
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Figure  13:  Receiver  Operating  Curves  for  the  DWBA  (dotted)  and  multiple  scattering  based 
(solid)  algorithms. 
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Figure  14:  Receiver  Operating  Curves  for  the  DWBA  (dotted)  and  multiple  scattering  based 
(solid)  algorithms. 
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Figure  15:  Receiver  Operating  Curves  for  the  DWBA  (dotted)  and  multiple  scattering  based 
(solid)  algorithms. 
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Figure  16:  Plots  of  Probabilities  of  detection  and  false  alarm  for  the  DWBA  (dotted)  and 
multiple  scattering  based  (solid)  algorithms  versus  threshold  value. 
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Figure  17:  Receiver  Operating  Curves  for  the  DWBA  (dotted)  and  multiple  scattering  based 
(solid)  algorithms. 
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Figure  18:  Receiver  Operating  Curves  for  the  DWBA  (dotted)  and  multiple  scattering  based 
(solid)  algorithms. 
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target  plus  clutter  K  matrix  is  measured  and  that  the  clutter  only  K  matrix  is  known  and 
the  detection  problem  then  reduces  to  making  a  decision  whether  or  not  target  is  present 
from  these  data.  The  signal  model  for  this  detection  problem  can  be  expressed  in  the  form 

(KtfM-KhM  +  WjtM  H\ 

Sj,k{Un)  =  <  f  (10) 

\W^k{un)  H0 

where  KSHC  is  the  experimentally  obtained  target  plus  clutter  K  matrix  at  5  Ghz,  Kfk  is  the 
experimentally  obtained  clutter  only  K  matrix  at  this  same  frequency  and  Wj^  is  AWGN 
that  is  generated  synthetically.  Here,  H\  is  the  hypothesis  that  target  is  present  and  Hq  is 
the  hypothesis  that  target  is  not  present. 


4.2  Detection  Schemes 


We  used  four  detection  schemes.  The  first  was  a  benchmark  that  made  its  decision  on  the 
total  signal  energy: 

f  choose  Hi  if  H^fcH  >T 

[choose  Hq  if  ||S,- *||  <T  1  ' 


where  T  is  a  threshold  that  is  varied  to  compute  the  ROC  curves  and  ||  •  ||  denotes  the 
Frobenius  norm  of  the  matrix.  The  benchmark  is  a  simple  energy  detection  scheme  and  does 
not  make  use  of  any  time  reversal  theory  of  methods. 

The  second  scheme  made  its  decision  based  on  the  maximum  singular  value  of  the  signal 
matrix;  i.e., 


choose  Hi  if  cq  >  T 
choose  Ho  if  cq  <  T 


(12) 


where  ai  is  the  first  (maximum)  singular  value  of  the  matrix  Sjtk-  This  scheme  is  in  the 
spirit  of  time  reversal  since  it  uses  the  SVD  of  the  measured  K  matrix. 

The  third  scheme  used  the  DORT  image  generated  from  the  singular  vector  corresponding 
to  the  largest  singular  value  of  the  Sjtk  matrix.  In  particular,  the  DORT  image  was  generated 
and  it  maximum  value  was  determined.  The  maximum  value  was  normalized  by  the  total 
energy  of  the  Green  function  vectors  at  the  image  point  and  this  statistic  was  compared  with 
a  detection  threshold  to  arrive  at  the  decision. 

The  final  scheme  was  similar  to  the  DORT  scheme  but  used  the  total  time  reversed  image 
field  that  is  generated  directly  from  the  measured  matrix  Sjtk-  This  image  field  was  discussed 
extensively  in  earlier  reports  and  briefings.  The  maximum  value  of  the  total  time  reversed 
image  field  was  normalized  by  the  total  energy  of  the  Green  function  vectors  at  the  image 
point  and  this  statistic  was  compared  with  a  detection  threshold  to  arrive  at  the  decision. 

It  is  important  to  note  that  the  last  two  schemes  require  that  the  location  of  the  target  be 
estimated  from  the  data.  The  results  of  the  study  indicate  that  this  is  very  important:  the 
SVD  of  the  K  matrix  is  not  enough.  It  is  important  to  use  a  test  statistic  that  is  tied  to  the 
estimated  value  of  the  time  reversed  image  fields  at  the  target  location. 
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Figure  19:  Real  data  simulation  #1  from  CMU  experimental  data. 
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Figure  20:  Real  data  simulation  #2  from  CMU  experimental  data. 
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Figure  21:  Real  data  simulation  #3  from  CMU  experimental  data. 
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5  Future  Research 

Whereas  the  research  performed  in  this  grant  and  reported  herein  concentrated  on  the  de¬ 
tection  of  point  targets  in  heavy  clutter  this  research  can  be  used  as  a  basis  for  extended 
target  imaging  and  inverse  scattering  applications.  In  particular,  the  generalized  MUSIC 
algorithms  developed  in  the  program  are  readily  extended  to  imaging  of  extended  targets 
while  the  inverse  scattering  algorithms  apply  with  virtually  no  change  to  these  more  general 
applications.  Candidate  applications  for  extension  of  the  results  obtained  in  the  project 
include: 

1.  Imaging  of  ground  targets  using  multistatic  data  collected  from  sets  of  UAV’s, 

2.  Imaging  in  urban  environments, 

3.  Imaging  of  hidden  and  unknown  targets  using  acoustic  data. 

The  above  research  goals  are  natural  extensions  of  much  of  what  has  been  performed  in 
the  program.  Moreover,  many  of  these  goals  were  originally  called  for  in  the  original  NU 
DARPA  proposal  and  much  of  the  groundwork  for  completing  these  goals  has  already  been 
laid  in  research  already  carried  out  in  the  project.  The  generalized  MUSIC  algorithms 
already  developed  in  the  project  can  be  cast  into  a  unified  conceptual  framework  based  on 
subspace  weighting  which  is  an  ideal  framework  for  achieving  goal  1.  This  formulation  leads 
naturally  to  a  tried  and  true  measure  or  performance  based  on  the  so-called  Point  Spread 
Function  (PSF)  of  the  imaging  algorithms.  Indeed,  this  measure  of  performance  and  the  use 
of  subspace  weighting  was  the  starting  point  of  the  famous  Backus  and  Gilbert  [14]  classical 
work  on  geophysical  data  inversion.  The  performance  of  the  inverse  scattering  algorithms 
can  also  be  based  on  computed  point  spread  functions  as  was  done  in  [8,  9,  10].  These  PSF’s 
are  functions  of  the  measurement  geometry,  bandwidth  and  antenna  radiation  patterns  and 
can  be  used  to  evaluate  system  parameters  in  terms  of  expected  performance. 

While  the  emphasis  in  goals  1  and  2  is  on  non-stochastic  or  “single  snap  shot”  data, 
in  Goal  3  the  emphasis  is  on  extending  the  theory  and  algorithms  developed  in  the  first 
two  goals  to  a  stochastic  environment.  Here,  the  performance  of  the  algorithms  would  be 
based  on  statistical  parameters  such  as  the  ensemble  average  PSF  and  the  use  of  Cramer 
Rao  performance  bounds.  Again  much  of  the  groundwork  for  this  goal  has  been  established 
in  the  project  [12,  13].  Further  extensions  would  include  the  development  of  statistically 
optimum  ML  estimation  algorithms  and  robust  estimation  and  imaging  algorithms  that  are 
immune  to  un-certainties  in  the  background  Green  function.  This  work  would  be  closely 
tied  to  background  Green  function  estimates  generated  in  Goal  2;  i.e.,  uncertainties  in  the 
computed  Green  function  will  be  employed  in  the  development  of  optimally  immune  imaging 
algorithms.  Some  preliminary  statistical  and  robustness  analysis  has  been  done  in  the  the¬ 
ses  of  Hanoch  Lev-Ari’s  students  Raha  Zandifar  (MS,  2001)  and  Minhtri  Ho  (PhD,  2005), 
including  peak  shape  for  TR-MUSIC,  effect  of  additive  noise  on  location  estimates,  and 
statistics  of  singular  values.  Much  more  remains  to  be  done,  especially  in  the  context  of 
non-ideal  wave  propagation  and  non-isotropic  antenna  elements. 
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